rm(list = ls())

## SET YOUR FILE PATH
setwd("~/Dropbox/Immigration/Analysis/Policy Diffusion/Replication Files/")

library(dplyr)
library(survminer)
library(tidyr)
library(standardize)
library(survival)
library(coxme)
library(frailtySurv)
library(frailtyEM)
library(frailtypack)
library(Hmisc)
library(haven)
library(foreign)
library(coefplot)
library(stargazer)
library(simPH)
library(coefplot)
library(xtable)
library(ggplot2)
library(texreg)
library(sp)
library(sf)
require(rgdal)
library(tmap)
library(tmaptools)
library(sp)
library(leaflet)
library(haven)
library(ggplot2)
library(shinyjs)
library(RColorBrewer)
library(ivtools)
library(jtools)

data <- read_dta("dwrap_index.dta")

data.africa <- subset(data, africa == 1)
data.onlytop <- subset(data, ihs_aidgdp_5ma_wins>=.0369172)
data.notop <- subset(data, ihs_aidgdp_5ma_wins<=.096386 )
data.1980 <- subset(data, year >= 1980)
data.1990 <- subset(data, year >= 1990)
data.indp5 <- subset(data, indp5 == 1)
data.indp10 <- subset(data, indp10 == 1)
data.law <- subset(data, numdomlaws >= 1)

set.seed(8675309)

## Set Convergence Criteria
coxph.control(eps = sqrt(.Machine$double.eps), toler.chol = .Machine$double.eps,
              iter.max = 1000, toler.inf = sqrt(.Machine$double.eps), outer.max = 1000)

coxme.control(eps = sqrt(.Machine$double.eps), toler.chol = .Machine$double.eps,
              iter.max = 1000)

## TABLE 4

res1 <- matrix(NA,35,9)

m1a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m1a)
res1[1:2,2] <- c(m1a$coef[1],sqrt(m1a$var[1,1]))
res1[4:5,2] <- c(m1a$coef[2],sqrt(m1a$var[2,2]))
res1[33,2] <- m1a$n
res1[34,2] <- m1a$history[[1]]$c.loglik
res1[35,2] <- AIC(m1a)

m2a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m2a)
res1[1:2,3] <- c(m2a$coef[1],sqrt(m2a$var[1,1]))
res1[4:5,3] <- c(m2a$coef[2],sqrt(m2a$var[2,2]))
res1[33,3] <- m2a$n
res1[34,3] <- m2a$history[[1]]$c.loglik
res1[35,3] <- AIC(m2a)

m3a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               daclib1last3+
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m3a)
res1[1:2,4] <- c(m3a$coef[1],sqrt(m3a$var[1,1]))
res1[4:5,4] <- c(m3a$coef[2],sqrt(m3a$var[2,2]))
res1[10:11,4] <- c(m3a$coef[11],sqrt(m3a$var[11,11]))
res1[33,4] <- m3a$n
res1[34,4] <- m3a$history[[1]]$c.loglik
res1[35,4] <- AIC(m3a)

keeps <- c("exc1inc2_reg_lag1", "std_ihs_aidgdp_5ma_wins", "std_ihs_fractal12_alesina", "landcontig", "ihs_gdppc_lag1", "ihs_pop_lag1", "std_polyarchy_lag1", "cwbroad_lag1", "ihs_iterate_lag1", "ihs_tradeshare_lag1", "majorannintra_contig_lag1", "numlib1", "fractalcontig", "govtfracprob", "lib1_start", "lib1_end", "liberalization1", "neggdpshock_lag1", "ccode2", "year")
df = data[keeps]

fitT.LX_aid <- coxph(formula=Surv(lib1_start, lib1_end, liberalization1)~
                       std_ihs_aidgdp_5ma_wins+
                       std_polyarchy_lag1+
                       exc1inc2_reg_lag1 +
                       majorannintra_contig_lag1+
                       ihs_gdppc_lag1 +
                       neggdpshock_lag1+
                       ihs_pop_lag1 +
                       cwbroad_lag1 +
                       ihs_iterate_lag1 +
                       ihs_tradeshare_lag1 +
                       
                       strata(numlib1)+
                       cluster(ccode2),
                     data = df,
                     method = c("efron"),
                     model = TRUE)
summary(fitT.LX_aid)

fitX.LZ_aid <- glm(formula=std_ihs_aidgdp_5ma_wins~
                     govtfracprob+
                     exc1inc2_reg_lag1+
                     std_polyarchy_lag1 +
                     ihs_gdppc_lag1 +
                     neggdpshock_lag1+
                     majorannintra_contig_lag1+
                     ihs_pop_lag1 +
                     cwbroad_lag1 +
                     ihs_iterate_lag1 +
                     ihs_tradeshare_lag1 +factor(ccode2), data=df)
summary(fitX.LZ_aid)

m4a <- ivcoxph(estmethod="ts", fitX.LZ=fitX.LZ_aid, fitT.LX=fitT.LX_aid, data=df, clusterid="ccode2", ctrl = F)
summary(m4a)
res1[1:2,5] <- c(-0.535, 0.061)
res1[4:5,5] <- c(-0.697, 0.059)
res1[33,5] <- c(2624)

m5a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins*std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m5a)
res1[1:2,6] <- c(m5a$coef[1],sqrt(m5a$var[1,1]))
res1[4:5,6] <- c(m5a$coef[2],sqrt(m5a$var[2,2]))
res1[7:8,6] <- c(m5a$coef[103],sqrt(m5a$var[103,103]))
res1[33,6] <- m5a$n
res1[34,6] <- m5a$history[[1]]$c.loglik
res1[35,6] <- AIC(m5a)

m6a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins*std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m6a)
res1[1:2,7] <- c(m6a$coef[1],sqrt(m6a$var[1,1]))
res1[4:5,7] <- c(m6a$coef[2],sqrt(m6a$var[2,2]))
res1[7:8,7] <- c(m6a$coef[170],sqrt(m6a$var[170,170]))
res1[33,7] <- m6a$n
res1[34,7] <- m6a$history[[1]]$c.loglik
res1[35,7] <- AIC(m6a)

m7a <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins*std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               daclib1last3+
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m7a)
res1[1:2,8] <- c(m7a$coef[1],sqrt(m7a$var[1,1]))
res1[4:5,8] <- c(m7a$coef[2],sqrt(m7a$var[2,2]))
res1[7:8,8] <- c(m7a$coef[104],sqrt(m7a$var[104,104]))
res1[10:11,8] <- c(m7a$coef[11],sqrt(m7a$var[11,11]))
res1[33,8] <- m7a$n
res1[34,8] <- m7a$history[[1]]$c.loglik
res1[35,8] <- AIC(m7a)

fitT.LX_aidint <- coxph(formula=Surv(lib1_start, lib1_end, liberalization1)~
                        std_ihs_aidgdp_5ma_wins*std_polyarchy_lag1+
                       majorannintra_contig_lag1+
                        exc1inc2_reg_lag1+
                       ihs_gdppc_lag1 +
                       neggdpshock_lag1+
                       ihs_pop_lag1 +
                       cwbroad_lag1 +
                       ihs_iterate_lag1 +
                       ihs_tradeshare_lag1 +
                       
                       strata(numlib1)+
                       cluster(ccode2),
                     data = df,
                     method = c("efron"),
                     model = TRUE)
summary(fitT.LX_aidint)

fitX.LZ_aidint <- glm(formula=std_ihs_aidgdp_5ma_wins~
                     govtfracprob+
                     std_polyarchy_lag1+
                     majorannintra_contig_lag1+
                       exc1inc2_reg_lag1+
                     ihs_gdppc_lag1 +
                     neggdpshock_lag1+
                     ihs_pop_lag1 +
                     cwbroad_lag1 +
                     ihs_iterate_lag1 +
                     ihs_tradeshare_lag1, data=df)
summary(fitX.LZ_aidint)

m8a <- ivcoxph(estmethod="ts", fitX.LZ=fitX.LZ_aidint, fitT.LX=fitT.LX_aidint, data=df, clusterid = "ccode2", ctrl = F)
summary(m8a)
res1[1:2,9] <- c(-2.455, 0.393)
res1[4:5,9] <- c(-0.388, 0.060)
res1[7:8,9] <- c(-0.525, 0.064)
res1[33,9] <- c(2624)

res1 <- round(res1,3)
res1[1,1] <- c("DAC Aid/GDP")
res1[4,1] <- c("Democracy")
res1[7,1] <- c("DAC Aid/GDP x Democracy")
res1[10,1] <- c("Top DAC Recipient Liberalization")
res1[33,1] <- c("Observations")
res1[34,1] <- c("Log-Likelihood")
res1[35,1] <- c("AIC")

latex(res1, file = "")
summary(m1a)
summary(m2a)
summary(m3a)
summary(m4a)
summary(m5a)
summary(m6a)
summary(m7a)
summary(m8a)

## FIGURE A.34

SimInt2 <- coxsimInteract(m6a, b1 = "std_ihs_aidgdp_5ma_wins", b2 = "std_polyarchy_lag1",
                          qi = "Marginal Effect",
                          X2 = seq(-.7, .6, by=.1),expMarg = F, nsim = 1000, ci = 0.95, spin=F, extremesDrop=T)

demintaid<-simGG(SimInt2, type = "points", alpha = .4, xlab = "\n Democracy (V-DEM)",
                 ylab = "Average Marginal Effect \n (AME) of DAC Aid/GDP \n", rug=T, lcolour = "#000000", pcolour = "#000000", method="loess")

demintaid<- demintaid + scale_x_continuous(breaks = seq(-.7, .7, by = .2)) + scale_y_continuous(breaks = seq(-1.5, 1, by = .25))

demintaid<-demintaid + geom_hline(yintercept=0, linetype="dashed", color = "black") + 
  ggtitle("Probability of Policy Liberalization")+
  theme(plot.title = element_text(hjust = 0.5))

demintaid

## Figure A.35: ADDITIONAL CONTROLS

m1b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               ustroopspc_lag1+
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
             frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data,
             method = c("efron"),
             model = TRUE)

summary(m1b)

afghanistan <- subset(data, country!="Afghanistan")
egypt <- subset(data, country!="Egypt")
iraq <- subset(data, country!="Iraq")
israel <- subset(data, country!="Israel")
jordan <- subset(data, country!="Jordan")
libya <- subset(data, country!="Libya")
pakistan <- subset(data, country!="Pakistan")

m2b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = afghanistan,
             method = c("efron"),
             model = TRUE)

summary(m2b)

m3b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = egypt,
             method = c("efron"),
             model = TRUE)

summary(m3b)

m4b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = iraq,
             method = c("efron"),
             model = TRUE)

summary(m4b)

m5b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = israel,
             method = c("efron"),
             model = TRUE)

summary(m5b)

m6b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = jordan,
             method = c("efron"),
             model = TRUE)

summary(m6b)

m7b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = libya,
             method = c("efron"),
             model = TRUE)

summary(m7b)

m8b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = pakistan,
             method = c("efron"),
             model = TRUE)

summary(m8b)

m9b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
               std_ihs_aidgdp_5ma_wins+
               std_polyarchy_lag1+
               exc1inc2_reg_lag1 +
               majorannintra_contig_lag1+
               ihs_gdppc_lag1 +
               neggdpshock_lag1+
               ihs_pop_lag1 +
               cwbroad_lag1 +
               ihs_iterate_lag1 +
               ihs_tradeshare_lag1 +
               
               cluster(ccode2) +
               strata(numlib1) +
               frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
               frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
             data = data.africa,
             method = c("efron"),
             model = TRUE)

summary(m9b)

m10b <- coxph(formula = Surv(lib1_start, lib1_end, liberalization1) ~
                std_ihs_aidgdp_5ma_wins+
                std_polyarchy_lag1+
                exc1inc2_reg_lag1 +
                majorannintra_contig_lag1+
                ihs_gdppc_lag1 +
                neggdpshock_lag1+
                ihs_pop_lag1 +
                cwbroad_lag1 +
                ihs_iterate_lag1 +
                ihs_tradeshare_lag1 +
                
                cluster(ccode2) +
                strata(numlib1) +
                frailty(ccode2, distribution = "gamma", method = "em", sparse = FALSE)+
                frailty(year, distribution = "gamma", method = "em", sparse = FALSE),
              data = data.onlytop,
              method = c("efron"),
              model = TRUE)

summary(m10b)

aidplot <- multiplot(m1b, m2b, m3b, m4b, m5b, m6b, m7b, m8b, m9b, m10b, secret.weapon=T,  innerCI=2, coefficients = c("std_ihs_aidgdp_5ma_wins"), sort = c("magnitude"), names=c("Control 1: US Troops/Capita", "Subset: No Afghanistan", "Subset: No Egypt", "Subset: No Iraq", "Subset: No Israel", "Subset: No Jordan", "Subset: No Libya", "Subset: No Pakistan", "Subset: Africa Only", "Subset: Excluding Country-Years \n in Bottom Quartile of DAC Aid/GDP"), horizontal = F, color = c("black", "black" , "black" , "black" , "black"), zeroColor ="red", textAngle = 0, title= "", xlab = "Standardized Coefficient  \n", ylab="")
aidplot
